/* Copyright 2001,2002,2003 NAH6
 * All Rights Reserved
 *
 * Parts Copyright DoD, Parts Copyright Starium
 *
 */
/*LINTLIBRARY*/
/*PROTOLIB1*/
//#include <math.h>
#include "main.h"
#include "bwexp.h"

/**************************************************************************
*
* ROUTINE
*               BWExpand
*
* FUNCTION
*               Bandwidth expansion of LPC predictor coefficients
*
* SYNOPSIS
*               BWExpand(alpha, pc, pcexp)
*
*   formal
*                       data    I/O
*       name            type    type    function
*       -------------------------------------------------------------------
*       alpha           fxpt_16 i       Bandwidth expansion factor
*       pc              fxpt_16 i       predictor coefficients
*       pcexp           fxpt_16 o       expanded predictor coefficients
***************************************************************************
*
* DESCRIPTION
*
*       Subroutine to perform bandwidth modification by moving the poles
*       (or zeros) radially in the z plane.  If the bandwidth expansion
*       factor (alpha) is less than unity, the bandwidths are expanded by
*       shifting the poles (or zeros) toward the origin of the z plane.
*       The predictor coefficients are scaled directly according to:
*
*                             i-1
*               a' = a  alpha           where i = 1, . . . , order+1
*                i    i
*
*       Resulting in a bandwidth expansion of:
*
*               -(fs/pi)ln(alpha) Hz
*
*       (e.g., fs = 8 kHz, alpha = 0.994127 -> 15 Hz bandwidth expansion)
*
*       CELP's LPC predictor coefficient convention is:
*              p+1         -(i-1)
*       A(z) = SUM   a   z          where a  = +1.0
*              i=1    i                    1
*
***************************************************************************
*
* CALLED BY
*
*       FindAdaptResidual FindImpulseResponse PostFilter
*       FindLPCResidual UpdateEverything LPC_Analysis
*
**************************************************************************/
void BWExpand(
fxpt_16 alpha,                  /* 0.15 format */
fxpt_16 pc[ORDER+1],            /* 2.13 format */
fxpt_16 pcexp[ORDER+1])         /* 2.13 format */
{
        int i;
        fxpt_16 powtab[ORDER];

FXPT_PUSHSTATE("BWExpand", -1.0, -1.0);
        powtab[0] = alpha;
        for (i=1; i<ORDER; i++)
                powtab[i] = fxpt_mult16_round(alpha, powtab[i-1], 15);

        pcexp[0] = pc[0];
        for (i = 1; i <= ORDER; i++)
                /*pcexp[i] = pc[i]*(float)pow((double)(alpha),(double)(i));*/
                pcexp[i] = fxpt_mult16_round(pc[i], powtab[i-1], 15);
FXPT_POPSTATE();
}

void BWExpand32(
fxpt_32 alpha,                  /* 0.31 format */
fxpt_32 pc[ORDER+1],            /* 5.26 format */
fxpt_16 pcexp[ORDER+1])         /* 2.13 format */
{
        int i;
        fxpt_32 powtab[ORDER];

FXPT_PUSHSTATE("BWExpand", -1.0, -1.0);
        powtab[0] = alpha;
        for (i=1; i<ORDER; i++)
                powtab[i] = fxpt_mult64_fix(alpha, powtab[i-1], 31);

//DUMPARR(powtab,31,ORDER)

//        pcexp[0] = fxpt_shr32_fast( pc[0], 13 );
        pcexp[0] = fxpt_saturate16_round( pc[0], 13 );
//        pcexp[0] = pc[0];
        for (i = 1; i <= ORDER; i++)
                /*pcexp[i] = pc[i]*(float)pow((double)(alpha),(double)(i));*/
//                pcexp[i] = fxpt_mult16_round(pc[i], powtab[i-1], 15);
//                pcexp[i] = fxpt_mult64_fix(pc[i], powtab[i-1], 28+16);
                pcexp[i] = fxpt_saturate16( fxpt_mult64_fix(pc[i], powtab[i-1], 28+16 ));
FXPT_POPSTATE();
}
